DSO 522: In-Class Workbook 4

Instructions

Your task is to edit this R Markdown notebook to include your solutions, ensuring that it is well-organized and professional in appearance. Answers must be clear and concise. Start by downloading the Rmd file.

Collaboration Policy. You are allowed to discuss this workbook with your classmates and work in groups. Despite group discussions, each student must write and submit their own solutions.

Assistance Policy. You may ask for clarifications from the instructor and teaching assistant in class. Do not seek help on these workbooks outside of class, as they are intended to be completed during class time.

Submission Requirements. Solutions must be submitted on Brightspace as a PDF writeup. Use the ‘Knit to PDF’ feature in RStudio to prepare your PDF document. Ensure that your PDF document looks like the provided PDF version of the workbook, including all code used to obtain the results in proper code blocks. Make sure your PDF writeup includes your name in the title (replace my name with yours).

Grading Criteria. This workbook is worth 10 points. To achieve a grade of 10, your writeup must correctly answer all questions, be easy to understand, and be formatted correctly.

Work Timeline. You are expected to work on this assignment in class, but you may complete it at home within 24 hours. It is due on August 28, 2024, at 9:30 pm. Late submissions will incur a 50% deduction for any initial delay below 24 hours, and an additional 10% deduction for each additional day.

Remember, you will need to use the fpp3 package to complete this assignment, as shown in the lecture.

suppressMessages(library(fpp3))

Problem 1

Produce forecasts for the following series using whichever of NAIVE(y), SNAIVE(y) or RW(y ~ drift()) is more appropriate in each case:

• Australian Population (global_economy)
• Bricks (aus_production)
• NSW Lambs (aus_livestock)
• Household wealth (hh_budget)
• Australian takeaway food turnover (aus_retail)

Problem 1 - Solutions

Australian population

global_economy |>
  filter(Country == "Australia") |>
  autoplot(Population)

Data has trend and no seasonality. Random walk with drift model is appropriate.

global_economy |>
  filter(Country == "Australia") |>
  model(RW(Population ~ drift())) |>
  forecast(h = "10 years") |>
  autoplot(global_economy)

Australian clay brick production

aus_production |>
  filter(!is.na(Bricks)) |>
  autoplot(Bricks) +
  labs(title = "Clay brick production")

This data appears to have more seasonality than trend, so of the models available, seasonal naive is most appropriate.

aus_production |>
  filter(!is.na(Bricks)) |>
  model(SNAIVE(Bricks)) |>
  forecast(h = "5 years") |>
  autoplot(aus_production)

Warning: Removed 20 rows containing missing values or values outside the scale
range (`geom_line()`).

NSW Lambs

nsw_lambs <- aus_livestock |>
  filter(State == "New South Wales", Animal == "Lambs")
nsw_lambs |>
  autoplot(Count)

This data does not appear to have much trend, although there may be some seasonality.

nsw_lambs |> filter(Month >= yearmonth("2016 Jan")) |> gg_season(Count)

Therefore, of the models available, seasonal naive is most appropriate.

nsw_lambs |>
  model(SNAIVE(Count)) |>
  forecast(h = "5 years") |>
  autoplot(nsw_lambs)

Household wealth

hh_budget |>
  autoplot(Wealth)

Annual data with trend upwards, so we can use a random walk with drift.

hh_budget |>
  model(RW(Wealth ~ drift())) |>
  forecast(h = "5 years") |>
  autoplot(hh_budget)

Australian takeaway food turnover

takeaway <- aus_retail |>
  filter(Industry == "Takeaway food services") |>
  summarise(Turnover = sum(Turnover))
takeaway |> autoplot(Turnover)

This data has strong seasonality and strong trend, so we will use a seasonal naive model with drift.

takeaway |>
  model(SNAIVE(Turnover ~ drift())) |>
  forecast(h = "5 years") |>
  autoplot(takeaway)

This is actually not one of the four benchmark methods discussed in the book, but is sometimes a useful benchmark when there is strong seasonality and strong trend.

The corresponding equation is \[ \hat{y}_{T+h|T} = y_{T+h-m(k+1)} + \frac{h}{T-m}\sum_{t=m+1}^T(y_t - y_{t-m}), \] where \(m=12\) and \(k\) is the integer part of \((h-1)/m\) (i.e., the number of complete years in the forecast period prior to time \(T+h\)).

Problem 2

Use the Facebook stock price (data set gafa_stock) to do the following:

1. Produce a time plot of the series.
2. Produce forecasts using the drift method and plot them.
3. Use an appropriate (but different) benchmark function to forecast the same data set.

Hint: The data must be made regular before it can be modeled. You can make the data regular, using trading days as index, as follows:

fb_stock <- fb_stock |>
  mutate(trading_day = row_number()) |>
  update_tsibble(index = trading_day, regular = TRUE)

Problem 2 - Solutions

a

fb_stock <- gafa_stock |>
  filter(Symbol == "FB")
fb_stock |>
  autoplot(Close)

An upward trend is evident until mid-2018, after which the closing stock price drops.

b

The data must be made regular before it can be modelled. We will use trading days as our regular index.

fb_stock <- fb_stock |>
  mutate(trading_day = row_number()) |>
  update_tsibble(index = trading_day, regular = TRUE)

Time to model a random walk with drift.

fb_stock |>
  model(RW(Close ~ drift())) |>
  forecast(h = 100) |>
  autoplot(fb_stock)

c

The most appropriate benchmark method is the naive model. The mean forecast is terrible for this type of data, and the data is non-seasonal.

fb_stock |>
  model(NAIVE(Close)) |>
  forecast(h = 100) |>
  autoplot(fb_stock)

Problem 3

Apply a seasonal naive method to the quarterly Australian beer production data from 1992. Check if the residuals look like white noise, and plot the forecasts. What do you conclude?

Hint: The following code will help.

# Extract data of interest
recent_production <- aus_production |>
  filter(year(Quarter) >= 1992)

Problem 3 - Solutions

# Extract data of interest
recent_production <- aus_production |>
  filter(year(Quarter) >= 1992)

# Plot the data
recent_production |>
  autoplot()

Plot variable not specified, automatically selected `.vars = Beer`

# Define and estimate a model
fit <- recent_production |> model(SNAIVE(Beer))

# Look at the residuals
fit |> gg_tsresiduals()

Warning: Removed 4 rows containing missing values or values outside the scale
range (`geom_line()`).

Warning: Removed 4 rows containing missing values or values outside the scale
range (`geom_point()`).

Warning: Removed 4 rows containing non-finite outside the scale range
(`stat_bin()`).

• The residuals are not centered around 0 (typically being slightly below it), this is due to the model failing to capture the negative trend in the data.

• Peaks and troughs in residuals spaced roughly 4 observations apart are apparent leading to a negative spike at lag 4 in the ACF. So they do not resemble white noise. Lags 1 and 3 are also significant, however they are very close to the threshold and are of little concern.

• The distribution of the residuals does not appear very normal, however it is probably close enough for the accuracy of our intervals (it being not centred on 0 is more concerning).

# Look at some forecasts
fit |>
  forecast() |>
  autoplot(recent_production)

The forecasts look reasonable, although the intervals may be a bit wide. This is likely due to the slight trend not captured by the model (which subsequently violates the assumptions imposed on the residuals).

Problem 4

Produce forecasts (for the next 5 years) for each of the 7 series from the state of Victoria in aus_livestock using SNAIVE(). Plot the resulting forecasts including the historical data. Is this a reasonable benchmark for these series? Do you see any issues?

Problem 4 - Solutions

aus_livestock |>
  filter(State == "Victoria") |>
  model(SNAIVE(Count)) |>
  forecast(h = "5 years") |>
  autoplot(aus_livestock)

• Most point forecasts look reasonable from the seasonal naive method.
• Some series are more seasonal than others, and for the series with very weak seasonality it may be better to consider using a naive or drift method.
• The prediction intervals in some cases go below zero, so perhaps a log transformation would have been better for these series.

Some limitations of this method are evident when we look at residuals. For example, this shows strong autocorrelation.

aus_livestock |>
  filter(State == "Victoria", Animal=="Bulls, bullocks and steers") |>
  model(SNAIVE(Count)) |>
  gg_tsresiduals()

Warning: Removed 12 rows containing missing values or values outside the scale
range (`geom_line()`).

Warning: Removed 12 rows containing missing values or values outside the scale
range (`geom_point()`).

Warning: Removed 12 rows containing non-finite outside the scale range
(`stat_bin()`).

Problem 5

Download daily prices for TSLA and GOOGLE, for the current year, from Yahoo!Finance using the following code.

library(tidyquant)
library(scales)

# Download the data from Yahoo!Finance
tickers <- c("TSLA", "GOOGL")
first_date <- "2024-01-01"
last_date <- Sys.Date()
prices.raw <- tq_get(tickers, get = "stock.prices", from = first_date, to = last_date)

# Because stock prices are not observed every day, we first set up a new time index based on the trading days rather than calendar days.
prices <- prices.raw |>
  group_by(symbol) |>
  mutate(day = row_number()) |>
  ungroup() |>
  as_tsibble(index = day, key=symbol) |>
  update_tsibble(index = day, regular = TRUE)

1. Produce a time plot of each of the two series.
2. Produce forecasts (10 trading days ahead) using the naive and drift methods, and plot them.

Problem 5 - Solutions

library(tidyquant)
library(scales)

tickers <- c("TSLA", "GOOGL")
first_date <- "2024-01-01"
last_date <- Sys.Date()

download_stock_prices <- function(tickers, first_date, last_date) {
  prices.raw <- tq_get(tickers, get = "stock.prices", from = first_date, to = last_date)

  prices <- prices.raw |>
    group_by(symbol) |>
    mutate(day = row_number()) |>
    ungroup() |>
    as_tsibble(index = day, key=symbol) |>
    update_tsibble(index = day, regular = TRUE)

  return(prices)
}

prices <- download_stock_prices(tickers, first_date, last_date)

a

prices |> autoplot()

Plot variable not specified, automatically selected `.vars = open`

b

prices |> 
  model(Naive=NAIVE(close), Drift=NAIVE(close~drift())) |> 
  forecast(h = 10) |> 
  autoplot(prices)

c

prices <- download_stock_prices("MSFT", first_date, last_date)

prices |> 
  autoplot(adjusted) +
  labs(title="Adjusted closing price of MSFT stock in 2024")

d

prices.ret <- prices |>
  arrange(date) |>
  mutate(return = adjusted / lag(adjusted) - 1)

prices.ret |>
  autoplot(return) +
  labs(title="Daily net returns for MSFT stock in 2024")

Warning: Removed 1 row containing missing values or values outside the scale
range (`geom_line()`).

Problem 6

Download daily prices for a different (not used in Problem 5) stock market symbol of your choice from Yahoo!Finance with tq_get()from the tidyquant package. Then answer the following questions using this data set.

1. Plot the time series of the ticker’s adjusted closing prices.
2. Compute daily net returns defined as \(r_t = p_t / p_{t-1} - 1\), where \(p_t\) is the adjusted day \(t\) price. Hint: the function lag() is helpful. Produce a time plot of the daily net returns. Do you see any patterns?
3. Produce forecasts for the daily net returns (10 trading days ahead) using the most appropriate method. Plot the forecasts along with the corresponding bootstrap prediction intervals.

Problem 6 - Solutions

library(tidyquant)
library(scales)

tickers <- c("MSFT")
first_date <- "2024-01-01"
last_date <- Sys.Date()

download_stock_prices <- function(tickers, first_date, last_date) {
  prices.raw <- tq_get(tickers, get = "stock.prices", from = first_date, to = last_date)

  prices <- prices.raw |>
    group_by(symbol) |>
    mutate(day = row_number()) |>
    ungroup() |>
    as_tsibble(index = day, key=symbol) |>
    update_tsibble(index = day, regular = TRUE)

  return(prices)
}

prices <- download_stock_prices(tickers, first_date, last_date)

a

prices |> 
  autoplot(adjusted) +
  labs(title="Adjusted closing price of MSFT stock in 2024")

b

prices.ret <- prices |>
  arrange(date) |>
  mutate(return = adjusted / lag(adjusted) - 1)

prices.ret |>
  autoplot(return) +
  labs(title="Daily net returns for MSFT stock in 2024")

Warning: Removed 1 row containing missing values or values outside the scale
range (`geom_line()`).

There is no visible pattern. There is also no significant autocorrelation.

prices.ret |> ACF(return) |> autoplot()

c

prices.ret |>
  model(MEAN(return)) |>
  forecast(h = 10, bootstrap = TRUE, times = 1000) |>
  autoplot(prices.ret)

Warning: Removed 1 row containing missing values or values outside the scale
range (`geom_line()`).